Compact binary evolutions with the Z4c formulation 
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Numerical relativity simulations of compact binaries with the Z4c and BSSNOK formulations 
are compared. The Z4c formulation is advantageous in every case considered. In simulations of 
non- vacuum spacetimes the constraint violations due to truncation errors are between one and three 
orders of magnitude lower in the Z4c evolutions. Improvements are also found in the accuracy of the 
computed gravitational radiation. For equal-mass irrotational binary neutron star evolutions we find 
that the absolute errors in phase and amplitude of the waveforms can be up to a factor of four smaller. 
The quality of the Z4c numerical data is also demonstrated by a remarkably accurate computation 
of the ADM mass from surface integrals. For equal-mass non-spinning binary puncture black hole 
evolutions we find that the absolute errors in phase and amplitude of the waveforms can be up to a 
factor of two smaller. In the same evolutions we find that away from the punctures the Hamiltonian 
constraint violation is reduced by between one and two orders of magnitude. Furthermore, the 
utility of gravitational radiation controlling, constraint preserving boundary conditions for the Z4c 
formulation is demonstrated. The evolution of spacetimes containing a single compact object confirm 
earlier results in spherical symmetry. The boundary conditions avoid spurious and non-convergent 
effects present in high resolution runs with either formulation with a more naive boundary treatment. 
We conclude that Z4c is preferable to BSSNOK for the numerical solution of the 3+1 Einstein 
equations with the puncture gauge. 



I. INTRODUCTION 

This is the conclusion of a series of papers (Tj-ll] about 
the development of a formulation of general relativity 
(GR), called Z4c, that attempts to combine the strengths 
of two popular evolution systems for applications in free- 
evolution numerical relativity. Here we summarize the 
logical development of the formulation. 

In the BSSNOK formulation there is a zero- 

speed characteristic variable in the constraint subsys- 
tem, which can result in large Hamiltonian constraint 
violations in numerical applications; the removal of this 
mode is one of the key advantages of the generalized 
harmonic formulation [9Ml2l] over BSSNOK. The gener- 
alized harmonic formulation also possesses a constraint 
damping scheme (ljjj . which exponentially damps away 
small, high-frequency constraint violations at the con- 
tinuum level. Furthermore the trivial wave-like nature of 
the generalized harmonic subsystem allows for the conve- 
nient construction of constraint preserving boundary con- 
ditions (l^ - fTrij . On the other hand, the key advantages 
in the BSSNOK formulation are the choice of confor- 
mal variables, and the fact that the formulation does not 
come tied to a particular gauge, which allows for the se- 
lection of the moving puncture gauge [l7l - l22j . The com- 
bination of these two strengths allows for the evolution of 
black holes represented by coordinate singularities on the 
grid without severe numerical difficulties. For BSSNOK 
radiation controlling constraint preserving boundary con- 
ditions have been proposed [23, but to our knowledge 
have not been successfully used in numerical applications. 
With these considerations in mind a conformal decompo- 



sition of the Z4 formulation [24l - l3fjj was proposed in 
Because of the close relationship between the Z4 and gen- 
eralized harmonic formulations, this conformal decompo- 
sition inherits all of the strengths outlined in the previous 
discussion. In [l[ a set of spherically symmetric (ID) tests 
involving single black hole and neutron star spacetimes 
demonstrated that, indeed, the Z4c formulation guaran- 
tees the robustness of, and better constraint preservation 
than BSSNOK, especially for non-vacuum spacetimes. It 
was also found that the main advantage of Z4c in the 
bulk of the computational domain, namely the propagat- 
ing constraint violations, presents problems at the outer 
boundary if the constraints are not absorbed but rather 
reflected. 

A first attempt to tackle this problem was presented 
in in which high derivative order constraint preserv- 
ing boundary conditions based on those of pH ] were pro- 
posed for Z4c. Well-poscdncss of the Z4c constraint 
subsystem initial boundary value problem was demon- 
strated, and numerical results in explicit spherical sym- 
metry demonstrated the efficacy of the constraint pre- 
serving boundary conditions. In this work we use con- 
straint preserving boundary conditions, motivated by a 
forthcoming study [f|, in which well-posedness of high 
order boundary conditions that are constraint preserving 
and control incoming gravitational radiation |3lL l32j , is 
analyzed. 

In [3[, the performance of the Z4 constraint damping 
scheme of [13( applied to Z4c in black hole and neutron 
star spacetimes was studied in detail in spherical sym- 
metry. The constraint damping scheme is effective, as 
expected, in the non-linear system provided that the con- 
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straint violation is sufficiently small and resolved on the 
numerical grid; in the case of grid noise the combination 
of artificial dissipation and damping helps to suppress 
constraint violations. But it was found that care should 
be taken in the choice of damping parameters. Success in 
spherical tests does not necessarily guarantee that of 3D 
simulations. Preliminary simulations with Z4c and Som- 
merfeld BCs in 3D with pure box in box mesh refinement 
showed poor behavior at the outer boundaries in long 
evolutions. Therefore to try and bridge the gap between 
the evolutions in spherical symmetry and full applica- 
tions in astrophysical spacetimes, numerical stability of 
Z4c in 3D evolutions was studied in Q , where numerical 
stability of the linearized Z4c system coupled to the punc- 
ture gauge with a novel discretization was shown, and 
numerical evidence from the Apples-with- Apples [33l435l | 
tests was presented with both the standard and the novel 
discretization. The importance of algebraic constraint 
projection was highlighted, and limitations of the punc- 
ture gauge for applications in cosmology were observed. 

The first application of Z4c was the study of the end- 
state of a collapsing neutron star with the puncture 
gauge [36| in I D and 3D simulations of spherical config- 
urations. The puncture gauge handles collapsing matter 
without the need for (matter or metric) excision because, 
during the process, the shift condition pushes the spatial 
coordinates off of the matter region. This study was fol- 
lowed up with a similar discussion of dust collapse in [37| ■ 

A variation of the conformal decomposition, CCZ4, 
was proposed in (38|. The difference between Z4c and 
CCZ4 is that CCZ4 includes parametrized constraint ad- 
dition that, for some choice, correspond to the origi- 
nal four-covariant Z4 formulation. In the generalized 
harmonic formulation it is known that the inclusion of 
these non-principal terms can be problematic in some 
test cases [39(, but constraint growth can be mitigated 
by the use of constraint damping, which is the approach 
of [38| . For a particular choice of the constraint addition 
parameters not corresponding to the four-covariant Z4 
system, evolutions of binary black hole spacetimes were 
presented and shown to reduce Hamiltonian constraint 
violation, by a factor of around four or five (see Fig. 4 
of |11|), relative to such simulations for the BSSNOK 
formulation. On the other hand, in the Z4c system non- 
principal constraint addition that makes the equations 
of motion as close as possible to those of the BSSNOK 
formulation whilst still obtaining the desired PDE prop- 
erties in the constraint subsystem is chosen. It would be 
interesting to know in more generality how the addition 
of non-principal constraints affects their evolution. 

In this paper we present the first long-term 3D evolu- 
tions of black hole and neutron star binaries with the Z4c 
formulation. In section [TT| the Z4c equations of motion 
and boundary conditions are presented. We then dis- 
cuss, in section [TTTl changes to the BAM numerical code 
since [4p| |. In particular we describe the implementation 
of spherical patches for the wave zone [4l| and the radi- 
ation controlling constraint preserving conditions of Q. 



In sections ITVl and IVl we present our simulations of single 
and binary compact objects, respectively. Appendices lAl 
and [Bl contain, respectively, descriptions of the spherical 
patch implementation and evolution of Teukolsky waves, 
the latter of which we use for code validation. We con- 
clude in section PVTl 

We use units G = c = 1 throughout, unless otherwise 
stated. 



II. THE Z4C EQUATIONS OF MOTION 

In this section we summarize for completeness the Z4c 
equations of motion, constraints, and boundary condi- 
tions. The evolved quantities of the formulation are the 
conformal spatial metric 7^, the lapse a, the shift vec- 
tor (3 l , the conformal tracefree part of the extrinsic curva- 
ture Aij , the constraint O and, up to constraint addition, 
the trace of the extrinsic curvature K = K — 29. Finally 
we evolve the the conformal contracted Christoffel sym- 
bols T l , which are initially set according to F l = —dffi. 

Evolution equations. The equations of motion for the 
Z4c formulation arc 



dtX= 3X [a(K + 2Q)-D l (3 t 
d t % = -2 a Aij + f3 k d k % + 2% {l d n (3 k 

for the metric components, 

d t K = -D l D ia + a A l3 A ij + -(K + 20) 2 

+ Ana [S + p] + am (1 - k 2 ) © + ftdiK . 
dtAij = xl-D.Dja + a (R tJ -8ir S^f 
+ a [(k + 2e)A ij -2A k i A kj 

+ f3 k d k A l3 + 2 A k(l d 3) p k - I A l3 d k p k , 

for the extrinsic curvature components and 
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for the remaining variables. Here the intrinsic curvature 
associated with the ADM metric 7^ = x _1 7ij is written 
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and we employ the shorthand 
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The derivative operator is that compatible with the 
ADM metric. Numerical evolutions are performed with 
a particular flavor of the puncture gauge [13, EM HH HH 



d t a = -o?il l K + f3 l dia, 



d t F = (X i iJisT'-riF + Fd j p 
Constraints. The system is subject to constraints 
= 0, 2Z l = f 



(11) 
(12) 



H = R + A l3 A lJ -\{K + 29) 2 - IQnp = 



(r d ) 4 = o, (13) 
), 

(14) 



djA ij + Y l 3k A ]k - r ' ! 0AK + 26) 



§^(logx)- 87^ = 0, (15) 
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of which the latter two, the algebraic constraints, are 
explicitly imposed in numerical integration. 

Boundary conditions. When the spacetime manifold 
has a smooth boundary with spacclike unit (with respect 
to the ADM metric) normal s* we choose for the trace of 
the extrinsic curvature the boundary condition 

d t K = - a JJTi (dj< + $K)- d A d A a + fdik , (17) 

where = denotes equality only in the boundary, we use 
the shorthand d s = s l di and we use upper case Latin 
indices A, B, C to denote those that have been projected 
with the operator q l j = 8 l j — s % Sj. We can alterna- 
tively write this as a second order derivative boundary 
condition on the lapse. Note that in the expressions for 
the boundary conditions we never commute the unit nor- 
mal, projection operator q l j, or physical projection op- 
erator q( ps>kl ij — q k iq l j — \qijq kl through any derivative 
operator. For example, dtT s = SidtT 1 . So we have 

8 t v i = 8%v a +q i A 9 t v A , (18) 
dtSij = SiSjdtS ss + ^qijdtSqq 4- 2s^qj^ A dtS sA 

+ q^ AB l3 d t S T A F B, (19) 



for vectors v l and symmetric tensors Sij, respectively. 
We refer to the various components of the time deriva- 
tive under this 2 + 1 decomposition as the scalar, vector 
and tensor sectors in the obvious way. If we are given 
conformally flat initial data with j3 l Si = and a 2 lis is 
constant, as is always the case for the data evolved in 
this work, we choose the boundary conditions 



d t T s =r 1 p l d i T s 
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for the longitudinal part of the shift. We take 
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for constraint preservation in the scalar sector, where for 
example D l = when acting on scalars, so that the 

tilde denotes that the conformal metric was used in the 
contraction. In the vector sector we have 



d t f A 



d s Y A - d A f s 

1 X. 
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for the gauge, and 
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for the constraints, where we denote projected indices by 
upper case characters starting from the beginning of the 
alphabet. Here we denote R qq = q^Rij. Finally we take 



d t A 



TF 

AB 



D S A AB - D {A A 



-A s{A b B) {\n X ) 



1 _ 9 - - i TF 

- A AB D s {\n X ) + A 1 A A iB - - A AB (K + 29) 



X D A Dl F a + C p A^ B 



TF 



(25) 



4 



for the tensor sector. The boundary conditions in the 
scalar and vector sectors are designed to absorb outgoing 
constraint violations. The last of the conditions (|25l) arc 
equivalent to the requirement that "Jo = (see [3ll.l32l.l42j 
for more details) . This condition could also be used with 
the BSSNOK formulation, but for constraint preservation 
more work may be needed to adapt the other conditions. 
In the development of this work we have tried alternative 
conditions on K and T l with only small differences in the 
outcome. We do not claim that these gauge boundary 
conditions arc optimal. In our numerical experiments 
we sometimes also employ the more naively constructed 
Sommerfeld conditions 

d t K=- ^Ea(d s K+±K) +P i d i K i (26) 
3 t f s =-^V^a(d s f s + ±f s )+/W s , (27) 
d t f A =-JlIs~a(d s f A + ±f A )+f3 l d l f A , (28) 

d t e=-a(d s e + },e) + p l d. l e, (29) 

9 t Ay= - a(d s A t] + i Aij) + $ h d k A\i ■ (30) 

Note that since the trace constraint on A y - is constantly 
imposed, the last of these conditions ([30)) constitute only 
five boundary conditions. Since the Z4c formulation cou- 
pled to the puncture gauge has ten incoming character- 
istics in the weak field region, these conditions are not 
overdetermined, in contrast to the standard conditions 
used with BSSNOK with box-in-box mesh refinement. 
Regardless of the formulation, or whether Sommerfeld 
conditions are taken for every evolved field or just the 
subset K,T l ,®,Aij, they are not constraint preserving 
and do not control incoming gravitational radiation. 

III. NUMERICAL METHOD AND 
PARAMETERS 

In this section we describe the numerical technique em- 
ployed in this work. We use the BAM code IH- 
[451 ] , a Cartesian-based adaptive mesh refinement (AMR) 
infrastructure optimized for the evolution of BBH and 
BNS spacctimes in 3+1 GR. Vacuum spacctime evolu- 
tions have also been performed with the AMSS-NCKU 
code (4^ | . which employs the same methods but with an 
independent implementation. 

AMSS-NCKU and BAM basics. Before discussing the 
upgrades to the codes used in this work, we summarize 
the main points of the numerical methods used by AMSS- 
NCKU and BAM. The evolution algorithm is based on 
the mcthod-of-lines with explicit Runge-Kutta (RK) time 
integrators (in this work we employed fourth order RK for 
vacuum spacetimes and third order RK for non-vacuum 
spacctimes) and finite differences approximation of the 
spatial derivatives. The numerical domain is made of a 
hierarchy of cell-centered nested Cartesian grids (nested 
boxes centered on the punctures [HI, H?} ) . The hierarchy 
consists of L levels of refinement labeled by I — 0, L— 1. 
A refinement level consists of one or more Cartesian grids 



with constant grid spacing hi on level /. A refinement fac- 
tor of two is used such that hi = h /2 l . The grids are 
properly nested in that the coordinate extent of any grid 
at level I, I > 0, is completely covered by the grids at 
level 1—1. Some of the mesh refinement levels can be 
dynamically moved and adapted during the time evolu- 
tion according to the technique of "moving boxes" . The 
Berger-Oliger algorithm is employed for the time step- 
ping [48j], though only on the inner levels (44|. Inter- 
polation in Berger-Oliger time stepping is performed at 
second order. A Courant-Fricdrich-Lcwy factor of 0.25 is 
employed in all the runs. We refer the reader to [4(], [46| 
for more details. 

Numerical treatment of the field equations. The BSS- 
NOK and Z4c equations of motion are implemented 
numerically in the same way. Fields derivatives are 
approximated by centered finite difference expressions 
(fourth order in this work), except for the shift advection 
terms which are instead computed with lop-sided expres- 
sions l49l - [5ll | . Algebraic constraints are enforced af- 
ter every time step in BAM and after every Runge-Kutta 
substcp in AMSS-NCKU. The gauge parameters in our 
numerical simulations are fixed to //£ = 2/a, /is = l/a 2 
and r\ — 2/Madm, unless otherwise stated. Note that 
this is not a choice of parameters for which the calcu- 
lations of H are expected to guarantee well-posedness 
of the initial boundary value problem because we have 
not carefully taken [ig in such a way as to avoid either 
some generically distinct speeds in the system clashing, 
or sets of measure zero on which the evolution equations 
may be only weakly hyperbolic. In earlier studies (38j | 
the choice c > with /is ~ c/a 2 was not found to 
greatly affect the behavior of the scheme in applications. 
As highlighted in Appendix A of Q it is challenging to 
identify problems caused by weak hyperbolicity in appli- 
cations, even if the system is weakly hyperbolic every- 
where in space. If the degeneracy happens on sets of 
measure zero we therefore expect that in practical ap- 
plications it will be nearly impossible to identify as the 
cause of any concrete numerical problem, although in 
principle such degeneracy should of course be avoided. 
In the Z4c simulations presented in this work the con- 
straint damping scheme with the values K\ = 0.02 and 
Ko = is used. These values have been suggested in 
the detailed ID numerical analysis of Q. Preliminary 
exploratory runs in 3D indicated the combined use of 
artificial dissipation and constraint damping terms is im- 
portant (in some cases essential) to avoid instabilities at 
the interfaces between boxes and spherical patches, thus 
confirming some of the expectations from the ID runs 
(see also the discussion in Appendix Ia"]) . OpenMP sup- 
port in terms of a hybrid OpenMP /MPI implementation 
has been added inside BAM to key functions with high 
computational cost, such as evolution equations, interpo- 
lation and wave extraction, which improves the efficiency 
of memory management of the code. 

Hydrodynamics treatment (BAM only). The algo- 
rithm implemented for the general relativistic hydro- 
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dynamics (GRHD, referred to hereafter as "matter" 
for brevity) is a robust high-resolution-shock-capturing 
(HRSC) method [43| based on primitive reconstruction 
and the Local-Lax-Friedrichs (LLF) central scheme for 
the numerical fluxes, see e.g. [52] . Metric variables 
are interpolated in space by means of sixth order La- 
grangian polynomials, while matter ones by a fourth or- 
der weighted-essentially-non-oscillatory scheme. Primi- 
tive reconstruction is performed here with the 5th order 
WENO scheme of [53[, which has been found important 
for long term accuracy [HI, HH . 

Spherical patches for the wave zone. Both the AMSS- 
NCKU and BAM codes have been upgraded for this 
work. The Cartesian box-in-box mesh refinement 
has been extended with spherical patches ("cubed 
spheres") [4l|, [HI, [57| for the wave zone. They provide 
us with adapted coordinates for waves, and, as demon- 
strated in appendix |A1 improved accuracy in GW extrac- 
tion. The presence of a spherical outer boundary further- 
more allows a straightforward implementation of BCs (ei- 
ther Sommcrfeld or CP) due to the absence of corners. 
When the spherical patches are being used the Carte- 
sian moving boxes, as previously implemented, are em- 
ployed only in the strong-field region for the simulation 
of the binary orbital motion, or a collapsing star, while 
the propagation of gravitational and possibly also electro- 
magnetic waves distant from the source is simulated on 
cubed spheres [4l|, H3] (see also [58l - [60l | ) . The technique 
is based on the covering of the sphere by six patches, each 
patch having local coordinates that are then mapped to 
Cartesian ones in such a way to avoid pathologies as- 
sociated with standard spherical polar coordinates. As 
opposed to the box-in-box setup, spherical patches al- 
low constant radial resolution with linear scaling in the 
number of grid points, while the boxes result in effec- 
tively constant angular resolution as well. In practice, 
the I = Cartesian box is substituted with six patches 
overlapping with the Cartesian box at level I = 1 and 
among themselves. The resolution of hi is also the radial 
resolution employed in the patches. A grid configuration 
is specified by the number of: i) levels L, ii) grid points 
in each non-moving box per direction n, iii) grid points 
in each moving box per direction n mv , iv) the coarsest 
box per direction h\ = h, v) grid points in each patch, n r 
and tiq ^ , which are typically chosen as ng t ^ = n/2. More 
details on the implementation are given in Appendix [Al . 
Finally, GWs are extracted using the Newman-Pcnrose 
formalism, in particular by computing the i/>4 scalar on 
coordinate spheres in the wave zone (see Sec. Ill of [ioj). 
Mode decomposition of -04 is performed by projections 
onto spin weighted spherical harmonics and integration 
on the spheres with a Simpson algorithm. 

Boundary conditions. The radiation controlling, con- 
straint preserving boundary conditions (|17H25p are im- 
plemented according to the following simple recipe. In- 
side every Runge-Kutta substep Lagrange extrapolation, 
of sixth order in our experiments, is used to populate 
enough ghostzones in a neighborhood of the boundary, 



so that the same finite difference and dissipation oper- 
ators used in the bulk may be evaluated at the bound- 
ary. The metric components a, p l ,Xilij arc updated at 
the boundary with their standard equations of motion 
from the bulk, whereas the boundary conditions (|17t - 
I25|) are used in place of the standard equations of mo- 
tion for K,T l ,Q, Aij. Since the evolution system is not 
symmetric hyperbolic, at least inside a large class of 
symmetrizers [j], with the standard puncture gauge, we 
can not rely on a discrete energy method to guaran- 
tee numerical stability, even in the linear constant co- 
efficient approximation. We will see however that this 
implementation of the boundary conditions is numeri- 
cally well-behaved in the experiments we perform. Som- 
merfeld boundary conditions are implemented in a sim- 
ilar way; instead of replacing the equations of motion 
for A',f l ,e,iij by $T7M$ we choose (HHl and like- 
wise for BSSNOK, but without the O boundary condi- 
tion. 



IV. LONG-TERM EVOLUTION OF SINGLE 
COMPACT OBJECTS 

In this section we present 3D evolutions of single iso- 
lated compact objects and compare systematically BSS- 
NOK and Z4c runs. The main focus is on the behavior of 
the Hamiltonian constraint violation, due to truncation 
and artificial BCs errors. We assess long term stability of 
the Z4c formulation in 3D evolutions of puncture black 
holes (both non-spinning and rapidly spinning) and com- 
pact stars, and demonstrate an overall improvement in 
constraint preservation and in some instances accuracy 
of physical quantities (absolute numerical errors at finite 
resolution) when the Z4c formulation is employed. The 
spurious effect of Sommcrfeld BCs and the improvement 
obtained with the new BCs are discussed in detail. Our 
results for spherically symmetric spacetimes are inter- 
preted in view of previous results obtained in 0, by 
means of ID simulations. The use of Z4c docs not signif- 
icantly improve the computation of the gravitational ra- 
diation emitted by a rapidly spinning Bowcn-York punc- 
ture [oll |. In particular we find that in both cases we do 
not obtain a clear pointwise convergence of the waves at 
the resolutions used in these tests. This result is possibly 
related to a lack of resolution, or to the use of punctures 
in the black hole description, rather than to deficiencies 
in cither formulation. 



A. Non-spinning puncture 

Evolutions of a non-spinning puncture of mass M = 1 
with BSSNOK and Z4c and Sommerfeld and constraint 
preserving BCs are compared. We confirm qualitatively 
the results obtained with ID simulations [lUf. In partic- 
ular the new BCs reduce significantly the spurious con- 
straint violation incoming from the boundary in the case 
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FIG. 1: Constraint violations in space at t = 75M (up- 
per panel) and t = 1000M (lower panel) in an evolution of 
a Schwarzschild puncture. At t = 75M we see an incoming 
pulse of Hamiltonian constraint violation in both the BSS- 
NOK and Z4c evolutions. The fact that the violation prop- 
agates in the BSSNOK evolution test is not in contradiction 
with the PDE properties of the BSSNOK constraint subsys- 
tem, because the Hamiltonian constraint itself is not a zero- 
speed characteristic variable of the constraint subsystem. At 
this resolution, the incoming constraint violation with the Z4c 
constraint preserving boundary conditions is roughly three 
times smaller than that of Z4c with the Sommerfeld bound- 
ary condition, but the violations with the constraint preserv- 
ing boundary conditions converge away with resolution whilst 
those of the Sommerfeld conditions do not. In the lower panel 
we can see the effect that the zero-speed mode of the BSS- 
NOK constraint subsystem has on the Hamiltonian constraint 
violation at the outer boundary at late times. 



of Z4c. The largest single constraint violation however 
occurs at the puncture, where both formulations give 
similar results. 

Setup. We employ five levels of box-in-box mesh re- 
finement, and attach the shells at r ~ 21.5 M. Each 
box has n = 40, and the resolution of level I = 5 
is he, = 0.0625 M per direction. We choose ne } $ — 40 
angular and n r = 40 radial points in each spherical patch 
so that the outer boundary is located at r = 58.5 M. No 
symmetries are imposed. The puncture is evolved with 
a precollapsed lapse as discussed in |4(j ■ In real applica- 
tions the outer boundary is typically placed further out, 
perhaps at 500 M or 1000 M. Although this does not 
solve issues caused by the boundary in principle, and is 
not an efficient treatment, in practice it reduces some of 
the features we will encounter here. By design however, 
in these tests, we aim to see the behavior of the constraint 
preserving conditions. 

Constraint violation in the strong-field region. In each 
test we find that the Hamiltonian constraint violation in 



a neighborhood of the puncture is large, taking a maxi- 
mum value, at this resolution of ~ 10 3 at around t = 5 M, 
before being rapidly suppressed to ~ 3. These values 
should not be taken particularly seriously because of the 
finite regularity of the solution at the puncture. Since 
much of the physics is concentrated around the puncture 
large numerical error and therefore large constraint vio- 
lation are to be expected in this region. These violations 
converge away more slowly than others in the simula- 
tions. The Z4c evolutions do not reduce this violation. 
As observed in spherical symmetry (see Fig. 10 of Q) we 
see small amounts of the puncture constraint violation 
propagating out of the horizon with either formulation. 
Unsurprisingly, the rectangular mesh-refinement bound- 
aries in the 3D evolutions obscure the feature. The Z4c 
evolutions do not qualitatively affect this behavior, which 
seems to improve with resolution. 

Constraint violation at large radii. In Fig. [T] the 
Hamiltonian constraint violation on the spherical shells 
in space at times t = 75 M and t = 1000 M is plotted. 
We find that the incoming constraint violation of the Z4c 
Sommerfeld evolution is comparable to that of the BSS- 
NOK evolution, although at late times most of the con- 
straint violation caused by the outer boundary has been 
absorbed by the boundary; the problem is that the con- 
straint violation induced by the Sommerfeld boundary 
condition, with cither BSSNOK or Z4c, does not converge 
away with resolution. Note that comparing the BSSNOK 
Sommerfeld evolution with the Z4c CPBC data is not 
entirely fair, because there is every possibility that con- 
straint preserving boundary conditions for BSSNOK, see 
for example (23|, could also be implemented. In any case, 
it is evident that the Z4c constraint preserving bound- 
ary conditions arc helpful in reducing this violation. It 
is possible to reduce further the constraint violation by 
using a large constraint damping term, as for example 
in [381 ] . but in spherical symmetry we found that values 
of K\ > 0.1 M can adversely affect the dynamics of the 
evolutions with Z4c at late times. We are therefore wary 
of using large damping parameters, but make no claim 
about the effect of such damping terms on other confor- 
mal decompositions of Z4. 



B. Non-rotating star 

Our comparison of the new three dimensional numerics 
with earlier findings in spherical symmetry is completed 
by the evolution of a single stable neutron star. We find 
Z4c advantageous in reducing constraint violation that 
accumulates in grid regions occupied by the matter. In 
particular the L2 norm of the Hamiltonian constraint is 
found to be three orders of magnitude smaller for stan- 
dard resolutions. The spurious ringing effect due to Som- 
merfeld BCs pointed out in spherical symmetry is visible 
also in 3D simulations, but does not dominate the er- 
ror budget at typical resolutions. In our earlier study [l[ 
we also evolved an unstable single star in a so-called mi- 
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gration test. We suppress such an experiment here; our 
tests of strong field dynamics are instead performed with 
binary spacetimcs in section IVl 

Setup. We use three levels of box mesh refinement, 
and attach the spherical grids at r ~ 18 M. The lowest 
resolution runs use n = 48, the star is completely covered 
by the innermost grid level I = 2 with a resolution of 
h,2 = 0.36 M per direction. We choose ng^ = 48 angular 
and n r = 48 radial points in each spherical patch so that 
the outer boundary is located at r = 50 M. Additional 
runs at twice the resolution (n = 96) are performed. No 
symmetries are imposed. The star is described by a T = 2 
polytropic EOS with gravitational mass M = 1.4 Mq and 
radius 5.7 M, which is exactly the same model evolved 
in [l[ in bespoke spherical symmetry. The evolutions 
are characterized by oscillations triggered by truncation 
errors. 

Constraint violation and truncation error dynamics. 
The key findings of the earlier numerical experiments 
in spherical symmetry were: (i) in the BSSNOK re- 
gion containing matter, the Hamiltonian constraint in- 
side the star grows throughout the evolution, whereas in 
the Z4c evolution it does not; (ii) the BSSNOK simu- 
lations have larger oscillations (larger truncation errors) 
than Z4c; (iii) if Sommerfeld BCs and sufficient resolu- 
tion are employed, incoming modes from the boundary 
perturb the star amplifying unphysically the oscillations 
(see Fig. 10). In 3D simulations we find some simi- 
lar features. The growth of the Hamiltonian constraint 
is the dominant one, because of lack of resolution the 
oscillations arc not significantly affected, and the effect 
of Sommerfeld BCs is less pronounced but visible. The 
norm of the Hamiltonian constraint is plotted as a func- 
tion of time in the left panel of Fig. [2l By the end of the 
simulation the constraint violation in the Z4c evolution 
is about two orders of magnitude smaller than the BSS- 
NOK violation and interestingly smaller even than the 
initial violation. The behavior is similar to that shown 
in Fig. 5 of [l|. We will demonstrate the growth of the 
Hamiltonian constraint also pointwise on the grid in sec- 
tion |Vl in the case of binary spacetimes. Considering 
momentum constraint violations one finds, as in spheri- 
cal symmetry, that the differences are far less dramatic. 
The main difference in this case is that the momentum 
constraint violations with BSSNOK are more dynamical, 
and slightly larger, although this is probably just because 
of reflections from the Sommerfeld outer boundary con- 
dition, which propagate many times over the computa- 
tional domain. With cither BSSNOK or Z4c, the largest 
persistent violation in the momentum constraint occurs 
at the surface of the star. In the right panel of Fig. [5] 
we show the oscillations of the central density during the 
simulation time. It is not possible to distinguish signif- 
icant differences, probably because of the low resolution 
employed. The large effect seen in [l[ was with a resolu- 
tion ten times higher than here. 

The Sommerfeld boundary kick. In the left panel of 
Fig. [Done sees the effect of the Sommerfeld BCs with Z4c 



on the norm of the Hamiltonian constraint. As in the sin- 
gle puncture evolution, roughly when the outer boundary 
becomes causally connected to the central body, there is 
a large incoming pulse of constraint violation. This pulse 
perturbs the central object, but unlike in the spherical 
case (see Fig. 1 of @), the violation is not the dominat- 
ing effect on the dynamics. In the right panel of Fig. [2] 
we do not see a significant jump in the central density. 
A blow-up of the central density plot is however shown 
in Fig. [31 together with data from a shorter run at twice 
the resolution. The figure demonstrates that at approx- 
imately 60M the star is slightly perturbed by the Som- 
merfeld BCs. More importantly, the figure shows that 
the size of the perturbation is not converging away with 
resolution, whereas the amplitude of the oscillations does, 
so we expect this error to be dominant at higher resolu- 
tions. By contrast CPBCs are characterized by smaller 
reflections. The Hamiltonian constraint violations prop- 
agating out from the star appear to converge at approxi- 
mately second order, in line with our expectation for our 
hydrodynamics scheme. This rate of convergence is main- 
tained by the constraint preserving boundary conditions 
(see also Fig. 5 of @). 



C. Rapidly spinning puncture 

Here we compare the evolutions of a single Bowen-York 
spinning puncture [6ll - l63l | with spin a = S/M 2 ~ 0.92. 
We choose such a comparatively large value of the spin 
parameter (sec also 64]) to test performance for large val- 
ues of the conformal factor. For a ~ 0.92 the puncture 
contributes only 24% of the mass, while the remainder 
is in the Brill wave contribution of the conformal factor. 
Both formulations give comparable results in terms of 
stability, the norms of the constraint violation, the ma- 
jority of which occurs at a few points near the puncture, 
and gravitational waves. With either system we observe 
that pointwise fourth order convergence is not achieved 
in the GWs at the resolutions at which we performed 
the tests, although global errors scale between third and 
fourth order. We expect that this behavior is related to 
inaccuracies caused either by lack of resolution, or intrin- 
sic to the puncture description of the black hole. 

Setup. We use five levels of box mesh refinement, and 
attach the spherical grids at r ~ 50 M. In the lowest res- 
olution runs each box has n = 48 points per direction, 
and the resolution of level I = 5 is h = 0.065 M. We 
choose = 48 angular and n r = 48 radial points in 
each spherical patch so that the outer boundary is lo- 
cated at r = 150 M. Runs at resolutions n = 64, 96, 128 
(with the grid spacing scaled in order to maintain the 
same grid setup) are performed. Bitant symmetry is im- 
posed. The initial data for this test is constructed in 
with an ad-hoc method, used elsewhere in the literature, 
in which the BAM spectral binary puncture initial data 
solver is applied to a single puncture with large spin and 
another, unboosted and unspinning, located very close 
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FIG. 2: The L2 norm of the Hamiltonian constraint violation (left) and the central density (right) as a function of time for 
evolution of a single stable star with polytropic EOS of index Y = 2. The Hamiltonian constraint violation is approximately two 
orders of magnitude smaller at the end of the experiment when using Z4c. With Z4c, the Sommerfeld boundary conditions (|26I 
I30|l cause large violation when the outer boundary becomes causally connected to the central body. The dynamics of all of the 
evolutions, the oscillations in the central density, are very similar. The star rings at its radial mode proper frequency. 




FIG. 3: A closer inspection of the oscillation of the central 
density at early times, both at the original and twice the 
resolution, for the Z4c evolutions. The constraint violation 
from the Sommerfeld boundary causes a jump in the central 
density, as observed in earlier work 1 , 2\ . This effect does not 
converge away with resolution, but at these resolutions is not 
the dominant effect. 



by with a relative uncorrected mass of 10~ 12 . There is 
no sign of the second puncture on the grid. Problems 
with convergence do however make this construction a 
point to address. 

Basic features of the dynamics. The (2, 0) multipo- 
lar mode of r ip 4 emitted during the evolution of this 



initial data are shown in Fig. @J for the lowest resolu- 
tion runs. The waves lie on top of each other; as reso- 
lution is increased they converge to each other. At late 
times, t = 350 M, a boundary effect is visible in the BSS- 
NOK data. The BCs of Z4c improves this behavior signif- 
icantly. Note that t = 350 M corresponds roughly to the 
time needed by a wave initially near the puncture to prop- 
agate to the outer boundary, be reflected from the strong 
field region and travel out once again to the extraction 
sphere at 90 M. No feature is visible as the wave passes 
from the outer boundary at around t = 210 A/, perhaps 
because ip 4 is by construction insensitive to incoming 
gravitational radiation. We computed the ADM mass in- 
tegral (see Eq. ([32} later) and find that in the BSSNOK 
runs a drift begins exactly when the outer boundary be- 
comes causally connected with each observer, although 
this drift is then swamped by the effect of the outgoing 
gravitational waves. The effect of the Sommerfeld bound- 
ary condition on physical quantities is discussed in more 
detail for binary neutron star simulations later. 

Constraint violation. As in the non-spinning case the 
largest Hamiltonian constraint violation is at the punc- 
ture. The evolution of the L2 norm of the Hamiltonian 
constraint is shown in the top panel of Fig. [5J for the low- 
est resolution runs. The results from Z4c and BSSNOK 
are comparable. At time t ~ 220 M BSSNOK data are 
affected by the Sommerfeld boundary conditions. The 
bottom panel of Fig. [5] shows the Hamiltonian violation 
in space at t ~ 220 M. The largest violation is at the 
puncture and of similar magnitude, but far from the grid 
origin Z4c violation behaves better. 

Convergence. We looked at self convergence of the 
waves presented in the upper panel of Fig. 0J At early 
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FIG. 4: Comparison of the waves with Z4c and BSSNOK 
from the single spinning puncture. The upper panel shows 
the (2,0) multipole mode of the real part of ri/>4 emitted 
during the evolution. The data is taken from the lowest reso- 
lution (48 3 points per direction) test. The Z4c and BSSNOK 
waves agree extremely well until about t ~ 300 M, roughly the 
time when the incoming constraint violation from the BSS- 
NOK Sommerfeld boundary condition is reflected from the 
central body and reaches the extraction sphere at r = 90 M. 
The lower two panels show the difference between the same 
quantity for the 64 (96) and 96 (128) point runs, respectively. 
The two systems are almost perfectly comparable. 
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FIG. 5: Hamiltonian constraint violation for the spinning 
puncture evolution. As in the upper panel of Fig. [4] data is 
taken from the lowest resolution (n — 48 points per direction) 
test. The upper panel shows the L2 norm of the constraint 
as a function of time. The jump in violation at around t = 
200 M in the BSSNOK data appears to be caused by the 
Sommerfeld boundary conditions, and does not converge away 
with resolution. The bottom panel shows the Hamiltonian 
violation in space at time t = 220 M, approximately at the 
peak of the jump in the upper plot, in the near-field (1 = 2 
box) region. 



times, up to t ~ 100, only second order pointwisc con- 
vergence is observed with either BSSNOK or Z4c. Later, 
the errors scale at third-to-fourth order rate in magnitude 
(norm) but pointwise convergence is lost. The use of ei- 
ther more resolution, properly constructed initial data, 
or simply more refinement levels, i.e. high resolution at 
the puncture, might improve this behavior. In (40| sim- 
ilar convergence tests were performed but with a lower 
spin a = 0.194, with two more, a total of seven, levels of 
mesh refinement, but with a smaller range in resolutions, 
and a lower maximum resolution in a neighborhood of 
the puncture. Although the plots presented in (see 
Fig. 6) are scaled for fourth order convergence, the differ- 
ence in resolutions makes it hard to distinguish between 
different convergence factors with confidence. The earlier 
study also found pointwise convergence is not maintained 
in the gravitational waves, thus our findings are consis- 
tent. A detailed discussion of spin and higher order finite 
differencing can be found in |64( . 



V. EVOLUTION OF COMPACT BINARY 
INITIAL DATA 

In this section we assess the performance of the Z4c for- 
mulation in the simulation of the merger of two compact 
objects. Two standard initial configurations are studied, 



with a set of approximately 60 simulations, in detail: an 
equal- mass, non spinning, black hole binary (BBH), and 
an equal-mass, irrotational neutron star binary (BNS). 
We discuss evolutions of about three orbits and compare 
systematically the results from several convergence series 
with the corresponding ones obtained with BSSNOK. 

Presentation rubric. In the presentation of the re- 
sults, we describe the main features of the dynamics, then 
the Hamiltonian constraint violations are compared and 
finally the accuracy of the GWs and other physical quan- 
tities are discussed in relation with the Hamiltonian con- 
straint violations. For brevity we discuss only the main 
emission channel, the (2, 2) multipole of the radiation. As 
is standard, the (complex) ip 4 projection ri/j% 2 , extracted 
on a coordinate sphere S r of radius r, is decomposed into 
amplitude and phase according to 

rV4=A 22 e-* 2 . (31) 
On the same coordinate sphere we compute the integral 

-EadmM = / ds / ^77 iJ 7 fe/ (-/ ik _j - 7y )ft ) , (32) 
Js r 

the spatial metric 7^ = X~ l lij which approximately rep- 
resents the ADM energy. The ADM energy (or mass) of 
the system is given by 

Madm = lim -Badm (?") , (33) 
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and it is a conserved quantity. On a given sphere in 
the wave zone, however, Eabm(t) is expected to devi- 
ate from Madm due to the gravitational energy radiated 
away from the sphere, 



-EradM = 



16tt 







j dt' J dn 


[ dt"ip 4 




Jo 



(34) 



Note that the outer integral in Eq. (f3"4")l is performed 
by a simple Riemann sum. The angular momentum is 
computed with similar ADM-like integrals, see e.g. [65^ 
and references therein, but note that this quantity re- 
mains ambiguously defined (gauge dependent) also in the 
asymptotic limit r — > oo for a generic asymptotically 
flat spacetime. We stress that the observed differences, 
up to outer boundary effects, seem to converge to the 
same continuum solution. They are however significant 
at fairly high resolutions in most of the cases we stud- 
ied. Particularly relevant arc the differences in Hamilto- 
nian constraint violation in the evolution of non- vacuum 
spacetimes. We show evidence that these violations are 
correlated with the quality of the numerical waveforms. 

Convergence tests. Standard three level self- 
convergence tests are performed using simulations at 
different grid resolutions. These tests can be biased 



by the choice of the triplet. As discussed in [55 1 
we consider only triplets in which i) the ratios be- 
tween the low and medium and medium and high 
resolutions arc h^/hnj — hm/hn > 1 (ideally > 2); 
and ii) the scaling factor is at least of order two, 
s = (h p L — h p M )/{h p M — h p H ) > 2, where p is the con- 
vergence rate. We find that if these criteria are not 
met then the measured convergence order from different 
triplets is not consistent. Additionally, in our experience 
we found it important to verify that the use of different 
triplets gives consistent results. Even when the criteria 
above are satisfied, different triplets can give differing 
convergence factors. The result is acceptable if the rates 
consistently improve as higher resolution triplets are 
taken. 

Initial data. Before discussing the results we sum- 
marize the initial data and grids employed. Initial 
data for BBHs are conformally flat puncture initial 
data constructed using the Brandt-Briigmann puncture 
method [66| with the BAM implementation of the spec- 
tral puncture initial data solver [671 ], The holes have 
equal mass, an initial separation of d = 7 M and are 
placed in a quasi-circular configuration, on which the 
Pade resummed eccentricity reduction algorithm of [68[ 
was applied. The initial data arc interpolated onto 
the AMSS-NCKU and BAM grids by eighth order La- 
grangian baryecntric interpolation. Initial data for BNS 
assume a conformally flat metric and irrotational flow. 
The initial separation is D ~ 10 M, ADM mass and an- 
gular momentum are M = 3.005 M Q and J = 8.3 M 2 , re- 
spectively. The stars are described by a T = 2 polytropic 
EOS, each has baryonic mass of Mf, = 1.625 Mq. These 
initial data have been produced with the LORENE |69[ li- 



TABLE I: Summary of the grid configurations and of the runs. 
Columns: name of the configuration, maximum refinement 
level, moving levels are those with I > i mv , number of points 
per direction in the moving levels, resolution in the level / = 
L— 1, number of points per direction in the non-moving levels, 
resolution in the levels 1 = 1 (radial resolution in the shells), 
number of radial points in the shells, number of angular points 
in the shells, outer boundary. Note that the resolution is 
given in units of M for BBH runs, but in units of Mq for 
BNS runs. Runs marked with "*" were reproduced with the 
AMSS-NCKU code. 



Name 


L 


L™ 




flL-l 


n 


hi 


n T 




Tout 


BBHO* 


9 


2 


48 


0.0182 


72 


2.33 


865 


24 


2092 


BBH1 


9 


2 


56 


0.0156 


84 


2.0 


1008 


28 


2091 


BBH2* 


9 


2 


64 


0.0137 


96 


1.75 


1150 


32 


2089 


BBH3 


9 


2 


72 


0.0122 


108 


1.556 


1293 


36 


2088 


BBH4* 


9 


2 


80 


0.0109 


120 


1.4 


1436 


40 


2088 


BBH5 


9 


2 


88 


0.0099 


132 


1.273 


1581 


44 


2090 


BBH6 


9 


2 


96 


0.0091 


144 


1.167 


1718 


48 


2087 


BBH7 


9 


2 


112 


0.0078 


168 


1.0 


2008 


56 


2088 


BBH8 


9 


2 


128 


0.0068 


192 


0.875 


2293 


64 


2086 


BBH9 


9 


2 


144 


0.0061 


216 


0.777 


2579 


72 


2086 


BNSO 


4 


1 


48 


0.5 


72 


2.0 


212 


24 


482 


BNSOr 


4 


1 


48 


0.5 


72 


2.0 


412 


24 


885 


BNS1 


4 


1 


56 


0.429 


84 


1.71 


245 


28 


482 


BNSlr 


4 


1 


56 


0.429 


84 


1.71 


478 


28 


882 


BNS2 


4 


1 


64 


0.375 


96 


1.5 


278 


32 


481 


BNS2r 


4 


1 


64 


0.375 


96 


1.5 


545 


32 


881 


BNS2a 


4 


1 


64 


0.375 


96 


1.5 


678 


32 


1081 


BNS3 


4 


1 


72 


0.333 


108 


1.33 


312 


36 


477 


BNS3r 


4 


1 


72 


0.333 


108 


1.33 


612 


36 


877 


BNS4 


4 


1 


80 


0.3 


120 


1.2 


345 


40 


476 


BNS4r 


4 


1 


80 


0.3 


120 


1.2 


678 


40 


875 


BNS5 


4 


1 


88 


0.273 


132 


1.09 


378 


44 


475 


BNS5r 


4 


1 


88 


0.273 


132 


1.09 


745 


44 


875 


BNS6 


4 


1 


96 


0.25 


144 


1.0 


412 


48 


476 


BNS6r 


4 


1 


96 


0.25 


144 


1.0 


812 


48 


876 



brary and have been already evolved in several places [43|, 
[55|,[ZQi. The data are interpolated onto the BAM grid by 
spectral interpolation. 

Grid setup. All of the runs performed at different res- 
olutions are listed in Tab. U together with details of the 
grid setup. Note that several of the BBH simulations, 
those marked with a "*" in Tab. HI have been performed 
with both the BAM and AMSS-NCKU codes, indicat- 
ing, at least for the chosen grid settings, the robustness 
of our findings. Note that our lowest resolution vacuum 
setup BBHO has maximum resolution comparable to that 
of the highest resolution in the earlier BAM calibration 
paper [40], although the boxes in [4(| were larger which 
also affects accuracy. The resolutions of the setups BBH2 
and BBH3 are comparable to the ones of recent BBH 
simulations [7lL [72[ that required only moderate accu- 
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racy, while the maximum BBH9 resolution is similar to 
what was used to obtain some of the BAM waveforms 
for the NINJA-2 catalog [z3|. The highest resolution in 
the matter simulations BNS6, BNS6r approaches those in 
the accurate runs of [HI, EH ■ Other details about gauge 
conditions, damping parameters, etc. have been already 
given in Sec. IIIII 



A. Equal-mass, non-spinning BBH 

Basic features of the dynamics. Let us discuss the 
evolution of BBH initial data. The black holes evolve 
for about 2.5 orbits before merging, radiating energy and 
angular momentum in gravitational waves. In Fig. [6] the 
puncture tracks are plotted for the grid BBH2 in Tab. fl] 
The gravitational waves have exactly the same basic pro- 
file, but unlike in the evolution of a single spinning punc- 
ture (see Fig. 2]) they are visually distinguishable, with 
the absolute maximum of the (2, 2) multipole mode of the 
real part of r ipi occurring, approximately 2.2 M, earlier 
in the BSSNOK data, in accordance with the expectation 
from Fig. [5J The differences that accumulate over the 
evolution converge away with resolution. In the BBH4 
data the delay is only 1.1 M. Also in contrast to the 
single spinning puncture case, presumably because the 
outer boundary is placed so far out at r ~ 2090 M, no 
boundary feature is visible in the BSSNOK waves, at 
least up to a radius of r = 400 M within the run-time of 
the simulation. 

Constraint violation. In Fig. [7] the Hamiltonian con- 
straint log 10 \H\ is plotted in space for the BBH1 runs 
on the orbital plane at a simulation time t — 146 M 
(roughly 1.5 orbits) on refinement levels I = 3,4,5. As 
in case of the single puncture evolution the Hamilto- 
nian constraint violation is dominated by the punctures 
and, for this data, at this time, has a maximum absolute 
value of ~ 10~ 2 regardless of the formulation. However, 
the Hamiltonian constraint violation differs in the strong 
field region depending on the formulation. In the BSS- 
NOK case a significant Hamiltonian violation extends to 
level I = 3 with an almost spherical pattern. In the Z4c 
case the violation is mainly restricted to the highest level 
around the puncture. Note the effect of Cartesian mesh 
refinement in the plot. For the same data the momentum 
constraint away from the puncture is also roughly an or- 
der of magnitude smaller in the Z4c data. The highest 
resolution runs, BBH9, have the smallest constraint vi- 
olation even in the strong field region, but the violation 
there is dominated by that of the puncture. The differ- 
ence between the Hamiltonian constraint violation of the 
BBH1 and the BBH9 data in this region is at most a fac- 
tor of three for either formulation despite the difference 
in resolution. 

Gravitational wave accuracy. Quantitative differ- 
ences are observed in the gravitational radiation com- 
puted in BSSNOK and Z4c simulations. The differences 
can be appreciated in self-convergence tests for phase 




-4 1 1 1 1 1 1 1 1 

-4-3-2-1 1 2 3 4 

x/M 

FIG. 6: Tracks of the punctures for binary black hole inspi- 
ral for the configuration BBH2 in Tab. [I] In the continuum 
limit, upto outer boundary effects, the tracks should agree 
perfectly because we are evolving the same data with the 
same gauge choice regardless of the formulation. At finite 
resolution however different formulations and discretizations 
may give different results. Initially the tracks agree, but a dif- 
ference accumulates; the BSSNOK punctures merge slightly 
sooner than the Z4c ones. This plot does not indicate either 
set of numerical data is better than the other, only that there 
is qualitative agreement between the results. 



and amplitude. Typical results are presented in Fig. [5] 
The phase errors accumulated to merger (t ~ 270 M) are 
shown in bottom panels. When the Z4c simulation is em- 
ployed we observe a factor three of improvement. This 
behavior is common to all of the simulations performed. 
Similarly, the amplitude error improves by a factor two- 
to-three, as can be appreciated from the top panel of the 
figure. The results obtained with the Z4c are not only 
characterized by smaller errors, but also by the fact that 
convergence is maintained for longer time, in particular 
in phase errors. All the errors are observed to converge 
with increasing resolution. The Z4c and BSSNOK data 
sets also appear to converge to each other, suggesting the 
same continuum limit is approached. We conclude that 
at finite resolution the use of Z4c results in the compu- 
tation of more accurate waveforms. 

Convergence issues. For the converge tests presented 
in Fig. [8] third order convergence is obtained, despite the 
use of fourth order operators for the bulk derivatives (see 
below). In our experiments with the described grid set- 
tings we find it difficult to achieve fourth order conver- 
gence in GWs from orbiting puncture runs, regardless 
of the formulation employed. On the other hand clearer 
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FIG. 7: Hamiltonian constraint violation in BBH simulations at time t ~ 146 M (after ~ 1.5 orbits) on the orbital plane and 
in the strong field region for the grid setup BBH1 in Tab. [T] The plots show log 10 \H\ for levels I — 5, 4, 3; the left panel shows 
the BSSNOK data, the right that of Z4c. Directly at the punctures the violation is similar in either case, but the surrounding 
region has smaller violation in the Z4c data. Some aspects of the grid structure are visible in the violation. 



(order approximately 3.5) convergence is demonstrated 
in Appendix [B] for the Teukolsky wave, a simpler prob- 
lem which however retains some of the important fea- 
tures, namely, nonlinearity, in the sense that we evolve 
with the full BAM Z4c infrastructure, constraint viola- 
tion, because the wave only satisfies the constraints to 
linear order in perturbation theory, present in the BBH 
calculations. On the other hand the Teukolsky wave test 
has neither moving boxes or punctures. We were un- 
able to identify the precise reason for the behavior in 
presence of punctures, however we point out the follow- 
ing well-known facts: i) punctures have finite regularity, 
which can obviously affect the formal order of conver- 
gence of high-order finite-differencing stencils; ii) moving 
box simulations have several source of errors (second or- 
der time interpolation, regridding, the number of mesh 
refinement buffer zones, 3D interpolation, data between 
boxes and shells, spheres and mode projection), a clearly 
defined error budget is very complicated to construct 
(for a discussion see e.g. [Io, EHzl])i hi) simulations 
have several freely specifiable parameters (grid parame- 
ters and/or gauge parameters), and systematically tun- 
ing all of them is beyond the aim of this work; iv) the con- 
vergence order observed actually depends on the triplets 
chosen in the self-converge test. The experimentally mea- 
sured convergence factors lie between two and four for 
both the formulations during the inspiral. Around the 
time of the merger the BSSNOK convergence rate drifts 
up to five or six. The drift is on average smaller for 
higher resolution triplets. This does not mean that the 
BSSNOK data are converging at a high order, but rather 
that the errors in the simulation do not allow us to judge 
the rate in a meaningful way with a simple error model. 



We do not claim that waveforms from either formulation 
are unreliable, but we are reluctant to estimate the abso- 
lute errors by Richardson extrapolation, which assumes 
a certain rate of convergence. 

AMSS-NCKU-BAM comparison. The qualitative 
features of these results have also been observed in 
simulations with the AMSS-NCKU code. Starting 
with the puncture tracks as in Fig. |5J we find that 
the BSSNOK also merge slightly earlier at a given 
resolution, roughly 3.2 M for BBH0 and 1.0 M for 
BBH2, in the same ball-park as the values obtained with 
BAM. Note that we do not expect to obtain identical 
values from the two codes, because although they share 
many ingredients, some specifics, for example the shells 
implementation and grid placement, differ. As in Fig. [7] 
the Hamiltonian constraint appears to be smaller in 
the Z4c data when using AMSS-NCKU. To test the 
robustness of this finding we have made a number of 
experiments with different precollapsed initial profiles 
for the lapse. We also tried different constraint damping 
factors for Ki, up to K\ = 0.1. Neither change seems 
to make a significant difference; the Z4c Hamiltonian 
constraint is always smaller. Note that in the earlier 
stability work |4[, the Hamiltonian constraint was 
computed directly from the conformal variables. Now, 
as in BAM, we compute it by transforming first to the 
ADM variables. The two versions of the calculation 
differ by additions of O and Zi, and obviously the 
finite difference approximation. The AMSS-NCKU 
code appears to produce slightly larger errors at mesh 
refinement boundaries than those of BAM. The reason 
for this is currently unclear, but may be the cause 
of the larger errors that we find in the AMSS-NCKU 
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FIG. 8: Convergence plot of binary black hole inspiral for the resolutions h — 1/56,1/80,1/112 (Runs BBH1, BBH4, and 
BBH7). The left panel shows results for BSSNOK, the right panel for Z4c. All the differences are scaled for 3rd order 
convergence. The extraction radius is at 150M. 



wave-forms. Nevertheless we still find, as in Fig. [5] that 
the Z4c wave forms are roughly twice as accurate in 
phase and amplitude. 



B. Equal-mass, irrotational BNS 

Basic features of the dynamics. Let us discuss the 
evolution of BNS. The binary evolves approximately two 
orbits before contact, then merges forming a hypermas- 
sive neutron star (HMNS) which finally collapses on dy- 
namical timescales. Dynamics and related gravitational 
wave emission were described in detail in our previous 
work [43| so we do not repeat them here. We only 
mention that the gravitational emission is characterized 
by approximately six cycles during which the GW fre- 



quency increases monotonically (even after contact), the 
merger time is defined conventionally at peak of the am- 
plitude's (2, 2) mode, nonlinear (quasi-radial and non- 
axisymmctric) oscillations of the HMNS generate the 
post-merger signal which decays exponentially after col- 
lapse. Figure [9] shows snapshots of the rest- mass density 
on the orbital plane around merger for a typical simula- 
tion obtained with the two formulation (and same setup) . 
There are visible differences in the "dynamics" , and in 
the Z4c run the contact and merger happens at a later 
simulation time. As in the case of the BBH puncture 
tracks, these plots are gauge dependent but they refer to 
the same gauge choice up to truncation errors and so they 
can be compared. In this respect, note that the centroids 
of the stars are offset, this is a coordinate effect due to 
the choice rj = 2 in the shift condition and was studied 
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FIG. 9: Evolutions of the rest-mass density on the orbital plane as computed with the grid setup BNS6 from Tab. [I] The results 
from the BSSNOK runs are plotted in the top panel and those from Z4c underneath. Similar comments to those in the caption 
of Fig. [5] apply; namely, it seems that the compact objects merge earlier in the BSSNOK data. In this case however there is 
already an expectation that the Z4c data will be more accurate because in earlier work [43| it was found that the objects merge 
later and later as resolution is increased with the BSSNOK formulation. 



in detail in [43|. 

Constraint violation. The L2 norm of the Hamilto- 
nian constraint violation during the evolution is reported 
in Fig. [TU] On the refinement level I = 2 one can observe 
an improvement of a factor ~ 100 during the whole evo- 
lution if Z4c is employed. By contrast the norms of each 
component of the momentum constraint agree for both 
evolutions, an almost constant violation around 10~ 7 is 
observed. The Hamiltonian violation on the orbital plane 
around merger time is shown in Fig. 1111 which includes 
now refinement levels I = 2 and 1 = 3. Even using the 
highest resolution run (BNS6) we register a two-to-three 
order magnitude difference in the absolute value of the 
Hamiltonian constraint. It is noteworthy that the viola- 
tion has an almost spherical pattern on and around the 
strong field region of the binary, which suggests that the 
violation is not dominated by rectangular mesh refine- 
ment boundaries. 

Gravitational wave accuracy. The differences at finite 
resolution between BSSNOK and Z4c described so far 
have an impact on the computation of physical quan- 
tities: gravitational waves and ADM mass. Figure [14] 
shows the result of a standard three-level self-convergence 
test for phase and amplitude of the (2, 2) mode of the 
gravitational radiation emitted. For the particular triplet 
shown BSSNOK data cease converging before contact. 
An analogue result was already found at similar resolu- 



tions in [43| (see also the detailed discussion in 0, [55[ ) . 
On the other hand the Z4c data are found to converge 
at approximately second order rate beyond contact and 
up to merger (t ~ 650 M). Similarly to the BBH case, 
Z4c waveforms arc found to be more accurate by a fac- 
tor three-to-four in accumulated phase and amplitude to 
merger time. We relate this behavior with the improve- 
ment obtained in Hamiltonian constraint preservation. 
Because truncation errors and the Hamiltonian violation 
(especially in the matter region) in BNS simulations are 
generically larger, the use of Z4c makes a significant dif- 
ference for the accuracy of these simulations. 

Mass and angular momentum conservation. As in the 
case of a single spinning puncture we consider the ADM 
mass integral, Eq. (p?2"j). and the energy radiated in GWs, 
Eq. (fM|). In Fig. [T3] we show that Z4c permits a reli- 
able computation of £admW- When corrected for the 
GW energy, the conservation of the ADM mass is of or- 
der of 0.1%. By contrast BSSNOK data do not even 
allow for a reliable estimate of Madm- Note that the 
ADM mass can be also estimated by means of a volume 
integral rather than a surface integral, see e.g. 
The more expensive volume integral computation can be 
more accurate and is found to give better results for BSS- 
NOK in the case of black hole binaries [7(| . Figure [T3] 
shows the ADM angular momentum integral, both with 
and without a correction by radiated angular momen- 
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FIG. 10: 1/2 Norm of the Hamiltonian constraint computed in 
the strong field, level I = 2, region, with the grid setup BNS6 
from Tab. U as in Fig. [5] At this resolution the Z4c violation is 
between one and two orders of magnitude smaller throughout 
the evolution. 



turn, at an extraction radius of 300 M, for both BSS- 
NOK and Z4c with CPBCs. The outer boundary for this 
run (BNS2a in Tab. IJ) is at approximately 360 M. Evi- 
dently at either radius the BSSNOK data is very poorly 
behaved, and has large error. In contrast in the Z4c 
evolution the corrected angular momentum are well con- 
served until the merger signal, which we expect to be 
much less accurate, reaches the extraction sphere. Early 
in the simulation the BSSNOK data are also conserved, 
but there is a jump at exactly the time t = 60 M when 
the outer boundary becomes causally connected to the 
extraction radius. We checked that this feature holds on 
every extraction sphere. This demonstrates clearly that 
the Sommerfcld boundary condition has an effect on the 
behavior of physical quantities inside BSSNOK simula- 
tions. One way to avoid the feature would be to place 
the outer boundary further away, which, depending on 
the physics of interest, may not be prohibitively expen- 
sive with the spherical shells for the wave zone, but it is 
not desirable to discard every extraction sphere as soon 
as it becomes causally connected to the outer bound- 
ary, because a larger domain would ideally be used for 
more reliable wave extraction, rather than as a buffer for 
poor boundary conditions. For comparison it would be 
very interesting to see data from BSSNOK simulations, 
where there is still a zero-speed mode in the constraints, 
with an implementation of the constraint preserving con- 
ditions of [231 ] . It is natural to compare the early times in 
Figures [T2l and ITUl to Fig. 11 of [73], a similar plot, com- 
puted with the spectral Einstein code, for inspiralling 
black hole neutron star binary data. The Z4c results are 
competitive, although it remains to be seen if they will 
continue to be so over many orbits. 



VI. CONCLUSIONS 

This paper is the conclusion of a scries [1 5] , the aim 
of which was to bring the advantages of the generalized 
harmonic formulation to the moving puncture method. 
We presented here for the first time 3D numerical rel- 
ativity simulations of compact binaries performed with 
Z4c, a conformal decomposition of the Z4 formulation. 

We started with evolutions of single compact objects 
and found that the expectations obtained by earlier work 
in spherical symmetry are largely borne out in the 3D nu- 
merics. The most striking feature in these tests is that, 
in the evolution of a single stable star by t = 1000 M at 
the resolutions used in our tests, is that the norm of the 
Hamiltonian constraint is approximately three orders of 
magnitude smaller in the Z4c data. For the first time 
we have presented results which combine radiation con- 
trolling, constraint preserving outer boundary conditions 
with the moving puncture method. These boundary con- 
ditions removed a perturbation to the central rest mass 
density of the star which is present when using Sommer- 
fcld boundary conditions. At the resolutions of our tests 
however the outer boundary is typically not the leading 
order contribution to numerical error. In evolutions of 
a single spinning puncture we found that the Z4c with 
the new outer boundary conditions remove certain con- 
straint violating features present in the BSSNOK data, 
but at least at early times the qualitative physical picture 
is unaltered. 

We then compared evolutions of compact binary space- 
times. In these tests we placed the outer boundary much 
farther away from the central body so as to simplify the 
discussion. Throughout the evolution of binary neutron 
star initial data, we find that at the same resolutions the 
Z4c formulation has between one and two orders of mag- 
nitude less Hamiltonian constraint violation in the norm. 
Interestingly the Hamiltonian constraint violation in the 
Z4c tests in this case stays at or below the level in the 
initial data, at least until the stars merge, and the simu- 
lations may therefore even be competitive with those of 
a constrained formulation. We find similar, albeit much 
less pronounced effects in the constraint violation in bi- 
nary black hole simulations, but the change of formula- 
tion does nothing to cure the dominant constraint vio- 
lation at the punctures themselves. The higher quality 
of the Z4c data is also apparent in physically meaning- 
ful quantities. In terms of gravitational wave accuracy 
we find that with any triplet, satisfying certain criteria, 
for either binary neutron star or binary black hole data, 
the absolute error in either amplitude or phase of the 
extracted gravitational waves is between two and four 
times smaller in the Z4c evolutions. The difference, in 
the evolution of compact binaries, between conservation 
of the ADM mass integral with the two formulations is 
remarkable. In the BSSNOK simulations one can not re- 
liably correct the integral with the radiated gravitational 
wave energy to arrive at a constant. In the Z4c simula- 
tions near perfect conservation is achieved. Furthermore, 
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FIG. 11: Hamiltonian constraint violation in BNS simulations at time t ~ 511 M (around merger) on the orbital plane and in 
the strong field region with the grid setup BNS6 from Tab. U as in Figs. I§1 and 1101 The plots show log 10 \H\; this plot is to the 
BNS data what Fig. is to the BBH data. In this case there is no sign of the Cartesian grid structure in the violation. Most 
of violation in the BSSNOK simulation appears inside the stars. The Z4c violation seems to be on average one or two orders 
of magnitude smaller than the BSSNOK violation, as seen more clearly in Fig. [9] 
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FIG. 12: ADM mass and radiated energy for BNS simulations with the grid setup BNS2a from Tab[I] Both are extracted inside 
the shells at R = 300M. The GW energy is computed from ip 4 according to [4(J]. The left panel shows results for BSSNOK, 
the right panel for Z4c; in comparison the latter demonstrates remarkable conservation. 



despite placing the outer boundary at a large coordinate 
radius, we find that the BSSNOK data are corrupted, 
for example in the angular momentum, by the Sommcr- 
feld boundary condition, whereas the Z4c data are free 
of this problem. If nothing else this motivates the use of 
constraint preserving boundary conditions with the BSS- 
NOK formulation. 



shown to give more accurate results than BSSNOK, both 
in terms of constraint violation and extracted physical 
quantities. We therefore expect that the Z4c formulation 
will become a standard tool for numerical relativity. 



In summary we have presented a large suite of nu- 
merical experiments in which the Z4c formulation was 
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FIG. 13: ADM angular momentum and radiated energy for BNS simulations with the grid setup BNS2a from Tab. [I] computed 
at an extraction radius of R = 300M. The radiated angular momentum Jaw is computed according to standard formulas, see 
for example [6^]. The left panel shows results for BSSNOK, the right panel for Z4c. The jumps in the BSSNOK data happen 
exactly as the extraction spheres become causally connected to the outer boundary at 360 M. Other extraction radii show the 
same feature. 
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spaced in local coordinates. The maps between the local 
coordinates and the Cartesian coordinates are 



±x patches 
±y patches 
±z patches 



arctan(y/x), 
arctan(a;/?/), 
arctan(a;/z), 



= arctan(z/a;) (Al) 
= arctan(z/y) (A2) 
= arctan(y/z) (A3) 



where (cj),9) e (— 7r/4 : 7r/4) x (— 7r/4 : 7r/4). The local 
radial coordinate R is a fisheye coordinate [78l - l80l | 



R = f(r) ~ /(0) , 



as function of the radius r = \/ x 2 + y 2 + z 2 . 
ticular function function f(r) implemented is 



(A4) 



The par- 



Appendix A: Spherical patches 

In this appendix the implementation of the spheri- 
cal patches in the BAM and AMSS-NCKU codes is de- 
scribed. We follow closely [U, to whom we refer for 
further details. We describe first the BAM implemen- 
tation and then highlight the differences in the AMSS- 
NCKU code. 

The grid structure. The BAM grid is made of a hi- 
erarchy of nested Cartesian grid boxes, each level is la- 
beled by the integer I = 0, 1, .... The level I = (outer- 
most box) is replaced with six patches with local coordi- 
nates Qj = {i?, 0, 9}. The grid of each patch is uniformly 



f(r) = Br + Aelog(cosh[(r-i?)/e)] (A5) 
d r f(r) = Atanh[(r-R)/e] +B , (A6) 

where the parameters A and B determine the stretching 
of the coordinate and e the transition region. The con- 
dition f(r) = r at interface between spherical patches 
and the boxes is important to avoid a step behavior and 
minimize interpolation errors. In this work we did not 
employ the fisheye coordinate, i.e. we set A = 0, B = 1. 

Spatial derivatives. Field derivatives at a grid point 
inside the patches arc calculated by finite differences on 
the uniformly spaced local grid. The derivatives in Carte- 
sian coordinates are obtained using the chain rule and the 
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FIG. 14: Convergence plot of binary neutron star inspiral for the resolutions h — 1/48,1/64,1/80 (runs BNS0, BNS2, and 
BNS4). The left panel shows results for BSSNOK, the right panel for Z4c. All the differences are scaled for 2nd order 
convergence. The extraction radius is at R = 400 M within the shells. The qualitative behavior of all curves does not change 
at different extraction radii. 



Jacobian of the transformation Xi(aj), 



d 

dxi 

d 2 



drj \ d 
dxi J drj ' 
drk dri \ d 2 



dxi dxj J drkdri \dxidxj J drk 



d 2 r k 







(A7) 



(A8) 



note that the last term in these equations fixes a typo in 
equation (5b) in j4lj . 

Grid schematic. The grid structure is sketched in 
Fig. [T5] The dark solid lines represent the physical grid, 
the lighter lines denote the ghost points which are needed 
for the finite differences. The green shaded regions de- 
note the ghost zones populated by inter-patch interpola- 
tion, and overlap with the neighbor patches. The other 



colored areas overlap with the / = level, and the res- 
olution of the Cartesian box is the same as the one of 
the radial one in the patch. The number of grid points 
in all the ghost zones (cyan, green, yellow) is equal. The 
distance between the points ro and r\ is equal to that 
between r\ and ri. 

Data communication. Spherical patches overlap on 
ghosts zones. Two neighboring patches share the radial 
coordinate and the angular coordinate perpendicular to 
the mutual boundary. Therefore only a ID interpolation 
parallel to the boundary has to be performed in the green 
regions of Fig. [T5] A Lagrange interpolation (sixth or- 
der in this work) which uses the most centered possible 
stencil is employed. Interpolation between patches and 
the I — level is performed with Lagrangian polynomial 
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FIG. 15: 2D scheme of the spherical patches. The darker 
regions denote the physical grids, the translucent numerical 
ghostzones. See the text for details of the various radii and 
how the color scheme relates to the numerical method. 



interpolation in 3D (colored regions of Fig. [T5|) . The var- 
ious interpolations are done in the following order: i) in- 
terpolate from box to shell (yellow region), ii) interpolate 
between shells (green region), hi) interpolate from shell 
to box (cyan and red regions), iv) set symmetries in box. 
For simplicity, grid symmetries are not applied in the 
shells during evolution. Each patch is evolved entirely 
and only afterward values at symmetry points are over- 
written by copying them. We note that because the in- 
terpolation of the red region in Fig. [TS] depends on points 
in the cyan region, the ordering of the different interpo- 
lation can give, in principle, different results. On the 
other hand, not interpolating the red region results in a 
"double evolution" , which we find leads to high frequency 
oscillations. 

Dissipation. Artificial dissipation is necessary for nu- 
merical stability during the evolution. In particular to 
maintain stable box-shell interface regions. Experimen- 
tally we found it important to apply different amounts of 
dissipation in box and spherical patches. In particular we 
use a lower dissipation on the spherical patches than in 
the box. For this work we tested only sixth order Kreiss- 
Oliger dissipation operators, using artificial dissipation 
coefficients a = 0.5 in the bulk and a = 0.1 in the shells 
(see [40| for our terminology). The number of angular 
points has been chosen according to the size of the I = 
box, ng tt p ~ n/2. Placing the box-shell interface close to 
the strong field region may also cause large growth of the 
error, even causing the code to crash, or affect the choice 
of the dissipation parameters required for stability. 

Parallelization. In order to reduce global MPI com- 
munication during shells synchronization, we use an opti- 



mized parallelization. The computations on the patches 
are distributed only in the radial direction. Every pro- 
cessor has six patch parts with the same radial exten- 
sion, hence the synchronization in the angular directions 
is performed locally by each MPI job. Note that this 
method implies a minimum number of radial points is 
necessary for a given number of processors. In some early 
tests we found that avoiding interpolation in this way can 
result in computations that are an order of magnitude 
faster in the shells. 

Variations of the numerical method in the AMSS- 
NCKU code. There are several differences in the im- 
plementation of shells in AMSS-NCKU with respect to 
BAM. In AMSS-NCKU there is no option to use a ra- 
dial fish-eye coordinate. Relating to the interpolation 
between box and shell, the scheme is somewhat different 
to BAM. We set the length between the point ro and r\ as 
the length of six points and we take these points as buffer 
points. We set the distance between r\ and r^ larger than 
the distance between r and r\. We take this part of re- 
gion as the double cover region of both box and shell. We 
set six points for the box outside of r2 and take them as 
buffer points. Similar to the treatment of the mesh refine- 
ment interface, we fill the buffer points for box and shell 
only at the end of a full RK4 step. As opposed to BAM 
we do not interpolate the points between r\ and r-i . Due 
to this double cover region, we can interpolate between 
box and shell in parallel fashion. But before that we have 
to synchronize the data between different shell patches. 
In the MPI communication, we divide the data in both 
radial direction and angular directions but try to make 
resulting blocks as cubic as possible for minimization of 
inter processor data change. 

Box-Shell comparison. Let us briefly discuss code 
performance and waveform quality by comparing runs 
which employ spherical patches ("shell runs" hereafter) 
in the wave zone to those that do not ( "box runs" here- 
after). For shell runs, radial and angular resolution can 
be adjusted separately, at linear scaling in the number of 
grid points in radius and quadratic scaling for both an- 
gular directions. Nested boxes imply effectively constant 
angular resolution for increasing radius and decreasing 
radial resolution. Increasing radial resolution leads to a 
cubic scaling in the number of grid points, because an- 
gular resolution increases simultaneously. This can be 
partially compensated by choosing larger boxes in the 
wave zone, which are comparatively cheap in run-time be- 
cause of Berger-Oliger time adaptivity, and typical runs 
are limited by run-time rather than memory. Constant 
radial resolution is helpful for tracking waves traveling to 
infinity, on the other hand decreasing radial resolution is 
a recipe sometimes employed deliberately for filtering fea- 
tures going to and coming from the outer boundary. An 
entirely different type of filter is effectively in use when 
extracting only (2, 2) modes. Features due to rectangular 
outer boundaries and rectangular refinement boundaries 
(cmp. Fig. [2]) are not visible at I = 2, but start to show 
at I = 4. A spherical outer boundary has the ad van- 
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FIG. 16: The BSSNOK box-shell comparison. The top row shows the errors in extrapolated waveforms when the wave zone is 
covered by Cartesian boxes as described in the text. The lower row shows the shell results. The left panel shows the amplitude 
and the right the phase of the (2, 2) mode of rifi 4 . The amplitude is improved by at least an order of magnitude and the phase 
by a factor around five. 



tage that a unique normal vector to the boundary exists, 
in fact our early (and incomplete) experiments with the 
new Z4c boundary conditions failed for a cubical outer 
boundary. However, a spherical boundary can also lead 
to a focused inward reflection of outgoing waves, while 
a rectangular boundary scatters spherical waves leading 
to a diffuse reflection. For example, the back-reflection 
feature in Fig. U is much smaller for cubical outer bound- 
aries, as are oscillations of a neutron star induced by 
boundary reflections. A clean treatment of a spherical 
boundary is nevertheless preferable since the reflections 
can be minimized, and even though a constraint violation 
from a cubical outer boundary may not be as visible, it 
is still present in the inner domain. The choice between 
box or shell runs actually depends on the waveform accu- 
racy goal, and requires a balance between accuracy and 
computational cost. We discuss two examples employing 
the BSSNOK formulation. 



BAM example. As a first example, we find that BAM 
BBH box runs with L = 9, L mv = 4, n mv = {48, 56, 64}, 
n = {82, 96, 110}, and h 9 = {0.0365, 0.0313, 0.0293}, 
give results comparable to shell runs with L — 7, 
L mv = 2i n mv = {48, 56, 64}, n = {82, 96, 110}, 
n r = {686, 800, 914}, n M = {14, 16, 18}, and h 6 = 
{0.0365, 0.0313, 0.0293}. Both series of simulations have 
the same resolution in the Cartesian boxes and a simi- 
lar computational cost. The resolution in the wave zone 
differs a factor two in the wave zone (lower in the box 
runs). The waveform errors and convergence properties 
are the same in both series. 

AMSS-NCKU example. As a second example we 
compare AMSS-NCKU BBH box runs with L = 
11, L mv = 8, n mv = {56,64 72}, n = {112, 128144} 
and h xl = {0.505/112, 0.505/128 0.505/144} M. With 
shell runs when the four coarsest levels are substituted 
by spherical patches that approximately extend to the 
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same outer radius. The computational cost of the shell 
runs is approximately twice the box runs. Using this 
triplet, third order self-convergence is achieved in the 
shell runs at all the extraction radii, while in the box 
runs self-convergence is lower than three at small radii 
and progressively degrades when the extraction is per- 
formed on coarser levels. The waveform in box runs is 
more noisy. Phase and amplitude of the waveform can 
be extrapolated according to 



K 



f(u,R) = J^f(u)R- k 



(A9) 



fc=0 



where f(u, R) is the quantity to extrapolate extracted at 
radius R, u is a retarded time, and R the Schwarzschild 
radius, see e.g. [HI, |8l|. A measure of the extrap- 
olation error is the difference between the extrapo- 
lated function and the last radius values. This er- 
ror is reported in the upper row of Fig. [TO] for box 
runs and the lower row for shell runs, using the radii 
R = {150, 140, 130, 120, 110, 100, 90, 80, 70, 60, 50} 
and different values of K, i.e. different polynomials. As 
evident from the figures the extrapolation for shell runs 
is more accurate, and results for different K are consis- 
tent with each other (no oscillations and overshooting). 
The "zig-zag" behavior in the left panel of the lower row 
of Fig. [TO] is due to the limited number of digits used in 
the output. For these data the choice K = 1 seems to be 
optimal. 



Appendix B: Evolution of Teukolsky waves 

In this appendix we evolve, for code validation, Teukol- 
sky's wave solution with the BAM code. 

Motivation. The codes used in this work have sources 
of error of different polynomial order in the grid spacing, 
but the operations performed most often, finite differenc- 
ing and time integration, are performed at fourth order. 
It is perhaps not surprising that we never find clear fourth 
order convergence in our evolutions of compact binaries. 

For the hydrodynamics simulations this behavior can 
be excused because the HRSC scheme we employ is only 
second order accurate. But for our vacuum simulations 
the situation is less clear. We therefore also studied the 
evolution of Teukolsky's solution [82| to the linearized 
Einstein equations. This solution represents a weak grav- 
itational wave propagating on a flat background. It sat- 
isfies the constraint equations of General Relativity at 
linear order. We use this solution as initial data and then 
evolve forward in time with the BAM code. Note that 
extensive convergence testing demonstrating fourth or- 
der convergence^in test cases, was shown for the AMSS- 
NCKU code in 



Teukolsky wave initial data. The metric has the form 

ds 2 = -dt 2 + (1 + Af rr )dr 2 + 2Bf r9 rdrd8 

+ 2Bf r4> rsm9drd<j) + (1 + Cfj$ + A,f (2 g >)r 2 d9 2 
+ 2{A - 2C)/ e0 r 2 smedOdc/) 

(Bl) 



+ (l + C/«+^y S in 2 ^ 2 



where we have 



and also 



f rr = sin #(cos </> — sin <f>) , 
f r e = sin0cos6'(cos 2 4> — sin 2 
f rcp — ~ 2 sin 9 sin <f> cos <p , 



/ e 7 = (l + cos 2 0)(cos 2 



A2) 

■lee 



f, 



(i) 



- (cos 2 4> - sin 2 4>) 

2 cos sin (f> cos </> , 

_M) 
Jtt ' 



fj?) = cos 9 cos 6>(cos 2 4> — sm ^ 



with 
A = 
B = 
C-- 



_3(^ 2 )^+3FW^ + 3F^ 



(>(3) ^ + 3f (2) ^ + 6F (1) ^ + gjF ^I 



4 V r 



2F 



(3); 



9F 



(2) 



21F {1) ^j + 21F^ 



(B2) 
(B3) 
(B4) 



(B5) 

(B6) 
(B7) 

(B8) 

(B9) 



(B10) 
, (Bll) 

(B12) 



This corresponds to an outgoing wave that is a pure m 
2 mode. The generating function F is given by 



(t-r) N , (t-r) 2 , 
F = « V - J cxp(-^-— 



X N 



X 2 



(B13) 



Derivatives of order n of F are denoted by F^ n K Apart 
from choosing a different generating function our wave is 
very similar to what was used in [83J. For our runs we 
have chosen a = 10 -4 , N = 10 and A = 10. 

Setup. We have evolved such a wave with an initial 
lapse of one and a shift of zero. For the evolution, as 
in the rest of the work, we have chosen the standard 
1+log and Gamma driver evolution equation for lapse 
and shift (|1HI12I) with fi L = 2/a,/j, s = I /a 2 . We use 
four levels of box mesh refinement, and attach the spher- 
ical grids at r ~ 30. In the lowest resolution runs each 
box has n — 48 points per direction, and the resolution 
of level I = 4 is h± = 0.167. We choose ng^ = n/2 an- 
gular and n r = n radial points in each spherical patch so 
that the outer boundary is located at r = 100. Runs at 
resolutions n = 48, 64, 80 (with the grid spacing scaled 
in order to maintain the same grid setup) are performed. 
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FIG. 17: Convergence plot for Teukolsky waves of the (2,2) mode of rip 4 scaled for fourth order convergence. Convergence 
series with n = {48, 64, 80}, and h 4 = {0.167,0.125,0.100}. The left panel shows the BSSNOK data, the right the Z4c. 
Extraction radius is at r — 80. 



Results. We find that both the BSSNOK and the Z4c 
system can successfully be used to evolve these waves. 
For Z4c we have set k,\ = 0.02 and «2 = 0. Since there are 
no strong gravitational fields we can extract the waves at 
any radius. Figure [17] shows a convergence plot for waves 
extracted at a radius of r = 80. The solid lines shows 
the difference between ifj 4 at low and medium resolution 
(L — M), while the dash dotted lines show the difference 



between medium and high resolution (M — H). When 
scaled with the proper factor of 3.66 for fourth order 
convergence, we see that M — H coincides approximately 
with L — M. For this data set the observed order of con- 
vergence is around 3.5. The left panel shows the results 
for BSSNOK, the right panel shows those for Z4c. 
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